#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "design.h"

// setup the value of the element [2*N+1], [2*N+2] of the array simplified
void setup_simplified (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k;

	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 0)
			continue;
        
		simplified[j][2 * N + 1] = simplified[j][2 * N + 2] = 0;
		for (k = n; k < N; k++)
		{
			simplified[j][2 * N + 1] += simplified[j][k];
			simplified[j][2 * N + 2] += simplified[j][N + k];
		}
	}
}

// apply the information of Fr to the simplified
void apply_Fr (int N, int L, int simplified[L][2 * N + 3], int Fr[N], int n)
{
	int j, k;

	/* Apply Forbidden R */
	for (k = 0; k < N; k++)
	{
		if (Fr[k] == 1)
		{
			for (j = 0; j < L; j++)
			{
				simplified[j][N + k] = 0;
			}
		}
	}

	/* Setup Array simplified */
	setup_simplified (N, L, simplified, n);
}

// apply the information Fr[p]=1 to the simplified
void apply_Fr2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int j;

	/* Apply Forbidden R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == -1 && simplified[j][N + p] == 1)
		{
			simplified[j][N + p] = 0;
			simplified[j][2 * N + 2] --;
		}
	}
}

// Find out all forbidden R
void simplify_Fr (int N, int L, int state[L][2 * N + 1], int simplified[L][2 * N + 3], int Fr[N], int i, int mode, int mode2)
{
	int j, k, flag;

	/* Get Forbidden R */
	for (j = 0; j < L; j++)
	{
		flag = 0;

		if (mode2 == 0)
		{
			flag = 1;
		}
		else if (mode2 == 1)
		{
			if (state[j][i] == state[j][N + i])
				flag = 1;
		}
		else if (mode2 == 2)
		{
			if (state[j][N + i] != state[j][2 * N])
				flag = 1;
		}
		else if (mode2 == 3)
		{
			if (state[j][i] == state[j][N + i] && state[j][N + i] != state[j][2 * N])
				flag = 1;
		}

		if (state[j][N + i] == 0)
		{
			simplified[j][2 * N] = -1;
			if (flag)
			{
				for (k = 0; k < N; k++)
				{
					if (state[j][k] == 1)
					{
						simplified[j][k] = 1;
						if (i != k)
							simplified[j][N + k] = 1;
					}
				}
			}
			else
			{
				for (k = 0; k < N; k++)
				{
					if (i == k)
					{
						if (state[j][k] == 1)
							simplified[j][k] = 1;
					}
					else
					{
						if (state[j][k] == 1 || state[j][N + k] == 1)
							simplified[j][k] = 1;
                                                
						if (state[j][k] == 1 && state[j][N + k] == 1)
							simplified[j][N + k] = 1;
					}
				}
			}
		}
		else
		{
			simplified[j][2 * N] = 1;
			if (flag)
			{
				for (k = 0; k < N; k++)
				{
					if (state[j][k] == 1)
					{
						simplified[j][k] = 1;
						if (i != k)
							Fr[k] = 1;
					}
				}
			}
			else
			{
				for (k = 0; k < N; k++)
				{
					if (i == k)
					{
						if (state[j][k] == 1)
							simplified[j][k] = 1;
					}
					else
					{
						if (state[j][k] == 1 && state[j][N + k] == 1)
							simplified[j][k] = 1;

						if (state[j][k] == 1 || state[j][N + k] == 1)
							Fr[k] = 1;
					}
				}
			}
		}
	}

	/* Apply Forbidden R and Reset Array simplified */
	apply_Fr (N, L, simplified, Fr, 0);
}

// apply the information of Fg to the simplified
void apply_Fg (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int Fg[N], int i, int n)
{
	int j, k;

	/* Apply Forbidden G */
	for (k = n; k < N; k++)
	{
		if ((i == k && Rr[k] == 1) || (i != k && Fg[k] == 1))
		{
			for (j = 0; j < L; j++)
			{
				if (simplified[j][k] == 1)
				{
					simplified[j][k] = 0;
					simplified[j][2 * N + 1] --;
					if (simplified[j][2 * N] == -1 && simplified[j][2 * N + 1] == 0)
					{
						simplified[j][2 * N] = 0;
					}
				}
			}
		}
	}
}
		
// apply the information Fg[p]=1 to the simplified
void apply_Fg2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int k;

	/* Apply Forbidden G */
	for (k = 0; k < L; k++)
	{
		if (simplified[k][2 * N] != 0 && simplified[k][p] == 1)
		{
			simplified[k][p] = 0;
			simplified[k][2 * N + 1] --;
			if (simplified[k][2 * N] == -1 && simplified[k][2 * N + 1] == 0)
				simplified[k][2 * N] = 0;
		}
	}
}

/* Return 1: find new member in Forbidden G */
int simplify_Fg (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int Fg[N], int i, int n)
{
	int j, k, flag;

	flag = 0;

	/* Get Forbidden G */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1 || simplified[j][2 * N + 2] != 0)
			continue;

		flag = 1;
		simplified[j][2 * N] = 0;
		for (k = n; k < N; k++)
		{
			if (simplified[j][k] == 1)
			{
				if (i == k)
					Rr[k] = 1;
				else
					Fg[k] = 1;
			}
		}
	}

	return flag;
}

// apply the information of Rg to the simplified
void apply_Rg (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Fr[N], int i, int n)
{
	int j, k, p;

	/* Apply Required G */
	for (k = n; k < N; k++)
	{
		if ((i == k && Fr[k] == 0) || (i != k && Rg[k] == 0))
			continue;

		for (j = 0; j < L; j++)
		{
			if (simplified[j][k] == 1)
			{
				if (simplified[j][2 * N] == 1)
				{
					simplified[j][2 * N] = 0;
				}
				else if (simplified[j][2 * N] == -1)
				{
					for (p = n; p < N; p++)
					{
						simplified[j][p] = 0;
					}
					simplified[j][2 * N + 1] = 0;
				}
			}
		}
	}
}

// apply the information Rg[p]=1 to the simplified
void apply_Rg2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int j, k;

	/* Apply Required R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != 0 && simplified[j][p] == 1)
		{
			if (simplified[j][2 * N] == 1)
			{
				simplified[j][2 * N] = 0;
			}
			else
			{
				simplified[j][2 * N + 1] = 0;
				for (k = 0; k < N; k++)
				{
					simplified[j][k] = 0;
				}
			}
		}
	}
}

/* Return 1: find new member in Required G */
int simplify_Rg (int N, int L, int simplified[L][2 * N + 3], int Rg[N], int Fr[N], int i, int n)
{
	int j, k, p, flag;

	flag = 0;

	/* Get Required G */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 1 && simplified[j][2 * N + 1] == 1)
		{
			flag = 1;
			simplified[j][2 * N] = 0;

			for (k = n; k < N; k++)
			{
				if (simplified[j][k] == 1)
				{
					if (i == k)
						Fr[k] = 1;
					else
						Rg[k] = 1;

					break;
				}
			}
		}
	}

	return flag;
}

// combination the R equations. For examples, the following two euqations:
//   r1r2(g3+g4)= 0 && r1r2(g5+g6)=0 can combine to r1r2(g3+g4+g5+g6)=0
void combinationR (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k, p, flag;

	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1)
			continue;

		for (k = j + 1; k < L; k++)
		{
			if (simplified[k][2 * N] != -1)
				continue;

			if (simplified[j][2 * N + 2] != simplified[k][2 * N + 2])
				continue;

			flag = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[j][N + p] != simplified[k][N + p])
				{
					flag = 1;
					break;
				}
			}

			if (flag == 1)
				continue;

			/* Starting combination */
			simplified[k][2 * N] = 0;

			if (simplified[j][2 * N + 1] == 0)
				continue;

			if (simplified[k][2 * N + 1] == 0)
			{
				simplified[j][2 * N + 1] = 0;

				for (p = n; p < N; p++)
					simplified[j][p] = 0;

				continue;
			}

			simplified[j][2 * N + 1] = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[k][p] == 1)
					simplified[j][p] = 1;
				simplified[j][2 * N + 1] += simplified[j][p];
			}
		}
	}
}

// simplify the R equation. For example, the following two equations:
// r1r2(g3+g4+g5+g6)=0 && g4+g5=1 can go to: r1r2=0 && g4+g5=1
void elimate1 (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k, p, flag;

	/* elimate R equation from G equation */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != 1)
			continue;

		for (k = 0; k < L; k++)
		{
			if (j == k)
				continue;

			if (simplified[k][2 * N] == 0)
				continue;

			if (simplified[k][2 * N + 1] < simplified[j][2 * N + 1])
				continue;

			flag = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[j][p] == 1 && simplified[k][p] == 0)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 0)
			{
				if (simplified[k][2 * N] == 1)
				{
					simplified[k][2 * N] = 0;
				}
				else
				{
					simplified[k][2 * N + 1] = 0;
					for (p = n; p < N; p++)
					{
						simplified[k][p] = 0;
					}
				}
			}
		}
	}
}

// apply the information of Rr to the simplified
void apply_Rr (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int i, int n)
{
	int j, k;

	/* Apply Required R */
	for (k = n; k < N; k++)
	{
		if (i == k || Rr[k] == 0)
			continue;

		for (j = 0; j < L; j++)
		{
			if (simplified[j][2 * N] == -1 && simplified[j][N + k] == 1)
			{
				simplified[j][2 * N] = 0;
			}
		}
	}
}

// apply the information Rr[p]=1 to the simplified
void apply_Rr2 (int N, int L, int simplified[L][2 * N + 3], int p)
{
	int j;

	/* Apply Required R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == -1 && simplified[j][N + p] == 1)
		{
			simplified[j][2 * N] = 0;
		}
	}
}

/* Return 1: find new member in Required R ; Return 2: find new member in Forbiddgen G ; Return 3: find both */
int simplify_Rr (int N, int L, int simplified[L][2 * N + 3], int Rr[N], int Fg[N], int n)
{
	int j, k, p, flag1, flag2;

	flag1 = flag2 = 0;

	/* Find Required R */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1)
			continue;

		if (simplified[j][2 * N + 2] == 1)
		{
			for (k = n; k < N; k++)
			{
				if (simplified[j][N + k] == 1)
					break;
			}

			p = k;

			if (simplified[j][2 * N + 1] == 0)
			{
				flag1 = 1;
				Rr[p] = 1;
				simplified[j][2 * N] = 0;
			}

			if (Fg[p] == 0)
			{
				flag2 = 1;
				Fg[p] = 1;

				/* Apply Forbidden G */
				apply_Fg2 (N, L, simplified, p);
			}
		}
	}

	if (flag1 == 1 && flag2 == 0)
		return 1;
	else if (flag1 == 0 && flag2 == 1)
		return 2;
	else if (flag1 == 1 && flag2 == 1)
		return 3;
	else
		return 0;
}

// simplify the R equation.  For example, the equation r1(g2+g3+g4)=0 imply
// the equation r1r3(g4+g4)=0.
void elimate2 (int N, int L, int simplified[L][2 * N + 3], int n)
{
	int j, k, p, flag;

	/* elimate R equation from R equation */
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] != -1)
			continue;

		for (k = 0; k < L; k++)
		{
			if (j == k)
				continue;

			if (simplified[k][2 * N] != -1)
				continue;

			if (simplified[j][2 * N + 2] > simplified[k][2 * N + 2])
				continue;

			if (simplified[j][2 * N + 1] == 0)
				;
			else if (simplified[k][2 * N + 1] == 0)
				continue;
			else if (simplified[j][2 * N + 1] < simplified[k][2 * N + 1])
				continue;

			flag = 0;
			for (p = n; p < N; p++)
			{
				if (simplified[j][2 * N + 1] != 0 && simplified[k][p] == 1 && simplified[j][p] == 0)
				{
					flag = 1;
					break;
				}

				if (simplified[k][N + p] == 0 && simplified[j][N + p] == 1)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 0)
			{
				simplified[k][2 * N] = 0;
			}
		}
	}
}

// set the array Og/Or based on the the array Rg/Rr/Fg/Fr, always return 0
int get_optional (int N, int L, int simplified[L][2 * N + 3], struct solve_static *solve1, struct solve_dynamic *solve2, int i, int n)
{
	int *Rg, *Rr, *Fg, *Fr, *Og, *Or;

	Og = solve1->Og;
	Or = solve1->Or;
	Rg = solve2->Rg;
	Rr = solve2->Rr;
	Fg = solve2->Fg;
	Fr = solve2->Fr;

	int j, k, flag;

	for (j = n; j < N; j++)
	{
		Og[j] = Or[j] = 0;
	}

	/* Get the option G */
	for (j = n; j < N; j++)
	{
		if (Rg[j] == 1 || Fg[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] == 0)
				continue;

			if (simplified[k][j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
			Og[j] = 1;
	}

	/* Get the option R */
	for (j = n; j < N; j++)
	{
		if (i == j || Rr[j] == 1 || Fr[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] != -1)
				continue;

			if (simplified[k][N + j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
			Or[j] = 1;
	}

	/* self inhibition */
	if (Rr[i] == 0 && Fr[i] == 0 && n <= i)
	{
		flag = 0;
		for (j = 0; j < L; j++)
		{
			if (simplified[j][2 * N] == 0)
				continue;

			if (simplified[j][i] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
			Or[i] = 1;
	}

	return 0;
}

// For irreducible solution only. Remove all removable interactions.
// Return value : 0 nothing is removed; 1 something is removed
int get_removable (int N, int L, int simplified[L][2 * N + 3], struct solve_dynamic *solve, int i, int n)
{
	int *Rg, *Rr, *Fg, *Fr;

	Rg = solve->Rg;
	Rr = solve->Rr;
	Fg = solve->Fg;
	Fr = solve->Fr;

	int j, k, flag, p;

	p = 0;

	/* Get the removable G */
	for (j = n; j < N; j++)
	{
		if (Rg[j] == 1 || Fg[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] == 1 && simplified[k][j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
		{
			p = Fg[j] = 1;
			apply_Fg2 (N, L, simplified, j);
		}
	}

	/* Get the removable R */
	for (j = n; j < N; j++)
	{
		if (i == j || Rr[j] == 1 || Fr[j] == 1)
			continue;

		flag = 0;
		for (k = 0; k < L; k++)
		{
			if (simplified[k][2 * N] == -1 && simplified[k][N + j] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
		{
			p = Fr[j] = 1;
			apply_Fr2 (N, L, simplified, j);
		}
	}

	/* self inhibition */
	if (Rr[i] == 0 && Fr[i] == 0 && n <= i)
	{
		flag = 0;
		for (j = 0; j < L; j++)
		{
			if (simplified[j][2 * N] == -1 && simplified[j][i] == 1)
			{
				flag = 1;
				break;
			}
		}

		if (flag == 0)
		{
			p = Fr[i] = 1;
			apply_Rg2 (N, L, simplified, i);
		}
	}

	return p;
}

// check the solution is irreducible or not.
// return value : 1 irreducible; 0 reducible
int check_irreducible (int N, int L, int simplified[L][2 * N + 3], struct solve_static *solve, struct solve_dynamic *solve2, int i)
{
	int j, k, flag, flag2, p;
	int Rg2[N], Rr2[N];
	int *Rg, *Rr, *edge;

	Rg = solve->Rg0;
	Rr = solve->Rr0;
	edge = solve2->edge;

	for (j = 0; j < N; j++)
	{
		Rg2[j] = Rr2[j] = 0;
	}

	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 0)
			continue;
		else if (simplified[j][2 * N] == 1)
		{
			if (simplified[j][i] == 1 && edge[i] == 0)
				continue;

			flag = 0;
			for (k = 0; k < N; k++)
			{
				if (k == i)
					continue;
				else if (simplified[j][k] == 1 && edge[k] == 1)
				{
					flag ++;
					if (flag == 1)
						p = k;
					else
						break;
				}
			}

			if (flag == 1)
				Rg2[p] = 1;
		}
		else if (simplified[j][2 * N] == -1)
		{
			if (simplified[j][2 * N + 1] > 0)
			{
				flag = 0;
				for (k = 0; k < N; k++)
				{
					if (k == i)
						continue;
					else if (simplified[j][k] == 1 && edge[k] == 1)
					{
						flag = 1;
						break;
					}
				}
                                
				if (flag == 0)
				{
					if (simplified[j][i] == 0)
						continue;
					else if (edge[i] == -1)
					{
						flag = 1;
						p = i;
					}
				}
				else
					flag = 0;
			}
			else
				flag = 0;

			for (k = 0; k < N; k++)
			{
				if (simplified[j][N + k] == 1 && edge[k] == -1)
				{
					flag ++;
					if (flag == 1)
						p = k;
					else
						break;
				}
			}

			if (flag == 1)
				Rr2[p] = 1;
		}
	}

	flag = 1;
	for (j = 0; j < N; j++)
	{
		if (edge[j] != 0 && Rg[j] + Rr[j] + Rg2[j] + Rr2[j] == 0)
		{
			flag = 0;
			break;
		}
	}

	return flag;
}

int check_solution (int N, int L, int simplified[L][2 * N + 3], int edge[N], int i, int old, int new)
{
	int j, k, ok, flag, flagR, flagG;

	flagR = flagG = 0;
	if (old == 1)
	{
		flagG = 1;
	}
	else if (old == -1)
	{
		flagR = 1;
	}
	else
	{
		if (new == 1)
			flagR = 1;
		else
			flagG = 1;
	}

	ok = 1;
	for (j = 0; j < L; j++)
	{
		if (simplified[j][2 * N] == 0)
			break;
		else if (flagG && simplified[j][2 * N] == 1)
		{
			flag = 0;
			for (k = 0; k < N; k++)
			{
				if (simplified[j][k] == 1 && k != i && edge[k] == 1)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 0 && (simplified[j][i] == 0 || edge[i] == -1))
			{
				ok = 0;
				break;
			}
		}
		else if (flagR && simplified[j][2 * N] == -1)
		{
			flag = 0;
			for (k = 0; k < N; k++)
			{
				if (simplified[j][N + k] == 1 && edge[k] == -1)
				{
					flag = 1;
					break;
				}
			}

			if (flag == 1)
				continue;

			if (simplified[j][2 * N + 1] == 0 || (simplified[j][i] == 1 && edge[i] == 0))
				flag = 1;
			else
			{
				for (k = 0; k < N; k++)
				{
					if (k == i)
						continue;
                                
					if (simplified[j][k] == 1 && edge[k] == 1)
					{
						flag = 1;
						break;
					}
				}
			}

			if (flag == 1)
			{
				ok = 0;
				break;
			}
		}
	}

	return ok;
}

// check the degree for the specified solution
int check_degree (struct designability *design, struct solve_dynamic *solve, int simplified[design->L - 1][2 * design->N + 3], int i)
{
	int N, L, *Rg, *Rr, *Fg, *Fr, *Og, *Or, *edge;

	N = design->N;
	L = design->L - 1;
	Rg = solve->Rg;
	Rr = solve->Rr;
	Fg = solve->Fg;
	Fr = solve->Fr;
	Og = solve->Og;
	Or = solve->Or;
	edge = solve->edge;

	int j, degree, tmp;
	int flag[3]; /* 0-R, 1-N, 2-G */

	degree = 0;

	for (j = 0; j < N; j++)
	{
		if (Rg[j] == 1 || Rr[j] == 1)
			continue;

		flag[0] = flag[1] = flag[2] = -1;

		flag[edge[j] + 1] = 1;

		if (Or[j] == 1)
			flag[0] = flag[1] = 1;
		else if (Fr[j] == 1)
			flag[0] = 0;

		if (Og[j] == 1)
			flag[1] = flag[2] = 1;
		else if (Fg[j] == 1)
			flag[2] = 0;

		if (design->dgflag != 2)
		{
			if (edge[j] == -1)
				flag[2] = 0;
			else if (edge[j] == 1)
				flag[0] = 0;
		}

		tmp = edge[j];

		if (flag[1] == -1)
		{
			edge[j] = 0;
			flag[1] = check_solution (N, L, simplified, edge, i, tmp, edge[j]);
		}

		if (flag[2] == -1)
		{
			if (flag[1] == 1)
			{
				edge[j] = 1;
				flag[2] = check_solution (N, L, simplified, edge, i, tmp, edge[j]);
			}
			else
			{
				flag[2] = 0;
			}
		}

		if (flag[0] == -1)
		{
			if (j == i && flag[1] == 1)
			{
				edge[j] = -1;
				flag[0] = check_solution (N, L, simplified, edge, i, tmp, edge[j]);
			}
			else
			{
				flag[0] = flag[1];
			}
		}

		edge[j] = tmp;

		flag[edge[j] + 1] = 0;
		degree += flag[0] + flag[1] + flag[2];
	}

	return degree;
}

// save or print or check degree for the specified solution
void get_solution (struct designability *design, struct solve_dynamic *solve, int simplified[design->L - 1][2 * design->N + 3], int i)
{
	int N, j, *edge;
	struct solution *solution;

	edge = solve->edge;
	solution = design->solution + i;

	N = design->N;

	if (design->svflag)
	{
		if (solution->max == solution->curr)
		{
			solution->max += 100;
			solution->edge = realloc (solution->edge, sizeof (int) * N * solution->max);
		}
        
		memcpy (solution->edge + solution->curr * N, edge, sizeof (int) * N);
		solution->curr ++;
	}
	else if (design->pnflag)
	{
		printf ("Node %2d :", i - design->mode1);
		for (j = design->mode1; j < N; j++)
		{
			if (design->mode1 == 1 && j == i && edge[0] == 1)
				printf (" g");
			else if (edge[j] == 0)
				printf (" n");
			else if (edge[j] == 1)
				printf (" g");
			else if (edge[j] == -1)
				printf (" r");
		}
		printf ("\n");
	}

	if (design->dgflag)
	{
		int degree;

		degree = check_degree (design, solve, simplified, i);
		design->degree_dist[i * (2 * N + 1) + degree] ++;
	}
}

// Enumerate all solutions and call get_solution for each solution 
void get_solution2 (struct designability *design, struct solve_dynamic *solve, int simplified[design->L - 1][2 * design->N + 3], int i, int n)
{
	if (n == design->N)
	{
		get_solution (design, solve, simplified, i);
	}
	else
	{
		if (solve->Og[n] == 1)
		{
			solve->edge[n] = 1;
			get_solution2 (design, solve, simplified, i, n + 1);
		}

		if (solve->Or[n] == 1)
		{
			solve->edge[n] = -1;
			get_solution2 (design, solve, simplified, i, n + 1);
		}

		if (solve->Og[n] == 1 || solve->Or[n] == 1)
		{
			solve->edge[n] = 0;
			get_solution2 (design, solve, simplified, i, n + 1);
		}
		else
		{
			get_solution2 (design, solve, simplified, i, n + 1);
		}
	}
}

void swap (int *a, int *b)
{
	int tmp;

	tmp = *a;
	*a = *b;
	*b = tmp;
}

int reorder (int N, int L, int simplified[L][2 * N + 3], struct solve_static *solve, struct solve_dynamic *solve2, int i, int j, int p)
{
	int *Rg, *Rr, *Fg, *Fr, *Og, *Or, *edge, *order;

	Rg = solve2->Rg;
	Rr = solve2->Rr;
	Fg = solve2->Fg;
	Fr = solve2->Fr;
	Og = solve->Og;
	Or = solve->Or;
	edge = solve->edge;
	order = solve2->order;

	int k;

	if (i == j)
		return p;
	else
	{
		for (k = 0; k < L; k++)
		{
			swap (&simplified[k][i], &simplified[k][j]);
			swap (&simplified[k][N + i], &simplified[k][N + j]);
		}

		swap (&Rg[i], &Rg[j]);
		swap (&Rr[i], &Rr[j]);
		swap (&Fg[i], &Fg[j]);
		swap (&Fr[i], &Fr[j]);
		swap (&Og[i], &Og[j]);
		swap (&Or[i], &Or[j]);
		swap (&edge[i], &edge[j]);
		swap (&order[i], &order[j]);
	}

	if (i == p)
		return j;
	else if (j == p)
		return i;
	else
		return p;
}

/************************************************************************
  Description of the input parameter:
  i : calculate the designability of the ith node
  n : enumerate the nth node (from 0 to n)
  flag : assign the value of flag to the nth node
************************************************************************/
int simplify (struct designability *design, struct solve_static *solve, struct solve_dynamic *solve1, int i, int n, int flag, int simplified[design->L - 1][2 * design->N + 3], int simplified0[design->L - 1][2 * design->N + 3])
{
	int N, L, j, k, p;
	int flag1, flag2, flag3;

	N = design->N;
	L = design->L - 1;

	struct solve_dynamic *solve2;
	int simplified2[L][2 * N + 3];

	solve2 = init_solve_dynamic (N);
	copy_solve_dynamic (N, solve1, solve2);

	memcpy (simplified2, simplified, sizeof (int) * L * (2 * N + 3));

	if (flag == 1)
	{
		solve->edge[n - 1] = 1;
		// set Required G flag to node (n-1)
		apply_Rg2 (N, L, simplified2, n - 1);
		// set Forbidden R flag to node (n-1)
		if (solve2->Fr[n - 1] == 0)
		{
			apply_Fr2 (N, L, simplified2, n - 1);
		}
	}
	else if (flag == 0)
	{
		solve->edge[n - 1] = 0;
		if (i == n - 1)
		{
			// set Required G flag to node (n-1) (no self-inhibition)
			apply_Rg2 (N, L, simplified2, n - 1);
		}
		else
		{
			// set Forbidden R flag to node (n-1)
			if (solve2->Fr[n - 1] == 0)
			{
				apply_Fr2 (N, L, simplified2, n - 1);
			}
			// set Forbidden G flag to node (n-1)
			if (solve2->Fg[n - 1] == 0)
			{
				apply_Fg2 (N, L, simplified2, n - 1);
			}
		}
	}
	else if (flag == -1)
	{
		solve->edge[n - 1] = -1;
		if (i == n - 1)
		{
			// set Forbidden G flag to node (n-1)
			apply_Fg2 (N, L, simplified2, n - 1);
		}
		else
		{
			// set Required R flag to node (n-1)
			apply_Rr2 (N, L, simplified2, n - 1);
			// set Forbidden G flag to node (n-1)
			if (solve2->Fg[n - 1] == 0)
			{
				apply_Fg2 (N, L, simplified2, n - 1);
			}
		}
	}

	// pseudo node's activation would forbid the self-inhibition
	if (design->mode1 == 1 && n == 1 && solve->edge[0] == 1 && solve2->Fr[i] == 0)
	{
		solve2->Fr[i] = 1;
		apply_Fr2 (N, L, simplified2, i);
	}

	// do the simplification
	if (flag != 2) // flag=2 means that we just setup an optional link, so no further simplification is required
	{
		flag3 = 0;
                
		do
		{
			p = 0;
			if (flag != -1)
			{
				p = simplify_Fg (N, L, simplified2, solve2->Rr, solve2->Fg, i, n);
				if (p == 1)
					apply_Fg (N, L, simplified2, solve2->Rr, solve2->Fg, i, n);
			}
                        
			flag1 = 0; // store the return value of get_removable
			flag2 = 2; // store the return value of simplify_Fr
                
			while (flag2 >= 2) // flag2>=2 means that we found new forbidden G, so we need to check required G and required R again
			{
				if (p == 1)
				{
					for (j = 0; j < L; j++)
					{
						if (simplified2[j][2 * N] == 1 && simplified2[j][2 * N + 1] == 0)
						{
							free_solve_dynamic (solve2);
							return 0;
						}
					}
				}
                
				if (flag != 1)
				{
					p = simplify_Rg (N, L, simplified2, solve2->Rg, solve2->Fr, i, n);
					if (p == 1)
						apply_Rg (N, L, simplified2, solve2->Rg, solve2->Fr, i, n);
				}
                
				if (flag != -1 || i == n - 1 || flag3 == 1)
				{
					/* Get Required R */
					if (flag1 == 0) /* Combine the R equation */
					{
						combinationR (N, L, simplified2, n);
					}
                                        
					elimate1 (N, L, simplified2, n);
                                        
					flag2 = simplify_Rr (N, L, simplified2, solve2->Rr, solve2->Fg, n);
					if (flag2 == 1 || flag2 ==3)
						apply_Rr (N, L, simplified2, solve2->Rr, i, n);
                                        
                                
					elimate2 (N, L, simplified2, n);
				}
				else
					break;
			}
                
			if (design->irflag == 0)
				flag1 = get_optional (N, L, simplified2, solve, solve2, i, n);
			else
				flag1 = get_removable (N, L, simplified2, solve2, i, n);

			flag3 = 1;
		}
		while (flag1); // flag1=1 means that some new links are marked as removable after calling the function get_removable, so need to check required G and required R again
	}

	if (n == 0)
	{
		memcpy (solve->Rg0, solve2->Rg, sizeof (int) * N);
		memcpy (solve->Rr0, solve2->Rr, sizeof (int) * N);
		for (j = k = 0; j < L; j++)
		{
			if (simplified2[j][2 * N] != 0)
			{
				memcpy (simplified0[k], simplified2[j], sizeof (int) * (2 * N + 3));
				k ++;
			}
		}
		for (j = k; j < L; j++)
			simplified0[j][2 * N] = 0;

		for (j = 0; j < N; j++)
		{
			if (solve->Og[j] == 1 && solve->Or[j] == 1)
				design->free_edges ++;
			else if (solve2->Rg[j] == 1)
				design->fix_edges ++;
			else if (solve2->Rr[j] == 1)
				design->fix_edges ++;
			else if (solve2->Fr[j] == 1 && solve2->Fg[j] == 1)
				design->fix_edges ++;
		}

	}

	if (n == N)
	{
		/* recover the order */
		for (j = 0; j < N; j++)
		{
			solve2->edge[solve2->order[j]] = solve->edge[j];
		}
        
		i = solve2->order[i];

		// if the solution is reducible, return 0
		if (design->irflag && check_irreducible (N, L, simplified0, solve, solve2, i) == 0)
		{
			free_solve_dynamic (solve2);
			return 0;
		}

		// calculate the total number of links and store to the variable p
		p = 0;

		for (j = 0; j < N; j++)
		{
			p += solve->edge[j] * solve->edge[j];
		}

		// update the information of minimum links
		if (p < design->min_link[i])
		{
			design->min_link[i] = p;
		}

		// setup the distribution of number of links
		if (design->bdflag && design->irflag)
			design->links_dist[i * (N + 1) + p] ++;
		else if (design->bdflag && !design->irflag)
		{
			int link_dist[N + 1], link_dist1[N + 1], link_dist2[N + 1], q, r;
                        
			link_dist[p] = 1;
			for (j = p + 1; j <= N; j++)
				link_dist[j] = 0;
                        
			for (j = 0; j < N; j++)
			{
				r = solve->Og[j] + solve->Or[j];
				if (r == 1)
				{
					memcpy (link_dist1, link_dist, sizeof (int) * (N + 1));
					link_dist2[p] = 0;
					for (k = p; k < N; k++)
					{
						link_dist2[k + 1] = link_dist1[k];
					}
					for (k = p; k <= N; k++)
					{
						link_dist[k] = link_dist1[k] + link_dist2[k];
					}
				}
				else if (r == 2)
				{
					memcpy (link_dist1, link_dist, sizeof (int) * (N + 1));
					link_dist2[p] = 0;
					for (k = p; k < N; k++)
					{
						link_dist2[k + 1] = 2 * link_dist1[k];
					}
					for (k = p; k <= N; k++)
					{
						link_dist[k] = link_dist1[k] + link_dist2[k];
					}
				}
			}
                        
			for (j = p; j <= N; j++)
				design->links_dist[i * (N + 1) + j] += link_dist[j];
		}

		// save or print the solutions
		if (design->svflag || design->pnflag || design->dgflag)
		{
			if (design->irflag)
			{
				get_solution (design, solve2, simplified0, i);
			}
			else
			{
				struct solve_dynamic *solve3; // save the information of array Rg/Rr/Fg/Fr/Og/Or/edge after the reodering

				solve3 = init_solve_dynamic (N);

				for (j = 0; j < N; j++)
				{
					solve3->Og[solve2->order[j]] = solve->Og[j];
					solve3->Or[solve2->order[j]] = solve->Or[j];
					solve3->edge[j] = solve2->edge[j];
				}

				if (design->dgflag)
				{
					for (j = 0; j < N; j++)
					{
						solve3->Rg[solve2->order[j]] = solve2->Rg[j];
						solve3->Rr[solve2->order[j]] = solve2->Rr[j];
						solve3->Fg[solve2->order[j]] = solve2->Fg[j];
						solve3->Fr[solve2->order[j]] = solve2->Fr[j];
					}
				}

				get_solution2 (design, solve3, simplified0, i, 0);

				free_solve_dynamic (solve3);
			}
		}

		free_solve_dynamic (solve2);
		return 1;
	}

	p = 1;

	int NS;

	if (design->mode1 == 1 && n == 0)
	{
		if (solve2->Rg[n] == 1)
		{
			solve->edge[n] = 1;
			NS = simplify (design, solve, solve2, i, n + 1, 2, simplified2, simplified0);
			free_solve_dynamic (solve2);
			return NS;
		}
	}
	else
	{
		for (j = n; j < N; j++)
		{
			if (solve2->Rr[j] == 1)
			{
				solve->edge[j] = -1;
				i = reorder (N, L, simplified2, solve, solve2, n, j, i);
				n ++;
			}
			else if (solve2->Rg[j] == 1)
			{
				solve->edge[j] = 1;
				i = reorder (N, L, simplified2, solve, solve2, n, j, i);
				n ++;
			}
			else if (solve2->Fg[j] == 1 && solve2->Fr[j] == 1)
			{
				solve->edge[j] = 0;
				i = reorder (N, L, simplified2, solve, solve2, n, j, i);
				n ++;
			}
			else if (design->irflag)
				;
			else if (solve->Og[j] == 1 && solve->Or[j] == 1)
			{
				solve->edge[j] = 0;
				i = reorder (N, L, simplified2, solve, solve2, n, j, i);
				n ++;
				p *= 3;
			}
			else if (solve->Og[j] == 1 && solve2->Fr[j] == 1)
			{
				solve->edge[j] = 0;
				i = reorder (N, L, simplified2, solve, solve2, n, j, i);
				n ++;
				p *= 2;
			}
			else if (solve->Or[j] == 1 && solve2->Fg[j] == 1)
			{
				solve->edge[j] = 0;
				i = reorder (N, L, simplified2, solve, solve2, n, j, i);
				n ++;
				p *= 2;
			}
		}
        
		if (n == N)
		{
			NS = p * simplify (design, solve, solve2, i, n, 2, simplified2, simplified0);
			free_solve_dynamic (solve2);
			return NS;
		}
		else if (n == N - 1)
			;
		else
		{
			int countG[N], countR[N], max, max_p;
                        
			for (j = n; j < N; j++)
			{
				countG[j] = countR[j] = 0;
				for (k = 0; k < L; k++)
				{
					if (simplified2[k][2 * N] == 0)
						continue;
					else if (simplified2[k][2 * N] == 1 && simplified2[k][j] == 1)
						countG[j] ++;
					else if (simplified2[k][2 * N] == -1 && simplified2[k][N + j] == 1)
						countR[j] ++;
                        
				}
			}
                        
			max = countR[n];
			max_p = n;
			for (j = n + 1; j < N; j++)
			{
				if (countR[j] > max)
				{
					max = countR[j];
					max_p = j;
				}
			}
			if (max > 0)
			{
				i = reorder (N, L, simplified2, solve, solve2, n, max_p, i);
			}
			else
			{
				max = countG[n];
				max_p = n;
				for (j = n + 1; j < N; j++)
				{
					if (countR[j] > max)
					{
						max = countG[j];
						max_p = j;
					}
				}
				i = reorder (N, L, simplified2, solve, solve2, n, max_p, i);
			}
		}
	}

	NS = 0;

	solve->Og[n] = solve->Or[n] = 0;

	if (solve2->Fg[n] == 0)
	{
		NS += simplify (design, solve, solve2, i, n + 1, 1, simplified2, simplified0);
	}

	NS += simplify (design, solve, solve2, i, n + 1, 0, simplified2, simplified0);

	if (solve2->Fr[n] == 0)
	{
		NS += simplify (design, solve, solve2, i, n + 1, -1, simplified2, simplified0);
	}

	free_solve_dynamic (solve2);
	return p * NS;
}

int check_satisfy (struct designability *design, int i)
{
	int j, k, N, L;

	N = design->N;
	L = design->L;

	int state[L - 1][2 * N + 1];

	/* setup the array state */
	for (j = 0; j < L - 1; j++)
	{
		for (k = 0; k < N; k++)
		{
			state[j][k] = design->pathway[j * N + k];
			state[j][k + N] = design->pathway[(j + 1) * N + k];
		}

		if (j != L - 2)
			state[j][2 * N] = design->pathway[(j + 2) * N + i];
		else
		{
			if (design->mode == 1)
				state[j][2 * N] = design->pathway[1 * N + i];
			else
				state[j][2 * N] = design->pathway[(j + 1) * N + i];
		}
	}

	/* special treatment for the trigger node */
	if (design->mode == 2)
		state[0][N + 2 * design->mode1] = 1;

	/* end of setup array state */

	int simplified[L - 1][2 * N + 3]; /* [0, N):g; [N, 2*N):r; 2*N:-1/0/1; 2*N+1/2:sum of g/r */
	int simplified0[L - 1][2 * N + 3];

	for (j = 0; j < L - 1; j++)
	{
		for (k = 0; k < 2 * N + 3; k++)
		{
			simplified[j][k] = 0;
		}
	}

	struct solve_static *solve;
	struct solve_dynamic *solve1;

	solve = init_solve_static (N);
	solve1 = init_solve_dynamic (N);

	/* No self-activation */
	solve1->Fg[i] = 1;

	/* Only activation is allowed for pseudo node(environment) */
	if (design->mode1 == 1) 
		solve1->Fr[0] = 1;

	/* Setup the links for the trigger node */
	if (design->mode == 2)
	{
		if (state[0][i] == state[0][N + i])        
		{
			solve1->Fr[design->mode1] = solve1->Fg[design->mode1] = 1;    // n     
		}
		else if (state[0][i] == 1)
		{
			solve1->Rr[design->mode1] = solve1->Fg[design->mode1] = 1;    // r
		}
		else
		{
			solve1->Fr[design->mode1] = solve1->Rg[design->mode1] = 1;    // g
		}
	}

	simplify_Fr (N, L - 1, state, simplified, solve1->Fr, i, design->mode, design->mode2);

	for (j = 0; j < L - 1; j++)
	{
		if (simplified[j][2 * N] == 1 && simplified[j][2 * N + 1] == 0)
		{
			free_solve_static (solve);
			free_solve_dynamic (solve1);
			return 0;
		}
	}

	design->NS2[i] = simplify (design, solve, solve1, i, 0, -2, simplified, simplified0);

	free_solve_static (solve);
	free_solve_dynamic (solve1);

	if (design->NS2[i] == 0)
		return 0;
	else
		return 1;
}

int calc_design (struct designability *design)
{
	int i, j, N;

	// Initialization
	N = design->N;
	design->NS = 1.0;
	design->min_links = 0;
	design->fix_edges = 0;
	design->free_edges = 0;

	if (design->bdflag)
	{
		for (j = 0; j < N * (N + 1); j++)
		{
			design->links_dist[j] = 0;
		}
	}

	if (design->dgflag)
	{
		for (j = 0; j < N * (2 * N + 1); j++)
		{
			design->degree_dist[j] = 0;
		}
	}

	if (design->svflag)
	{
		for (j = 0; j < N; j++)
		{
			design->solution[j].curr = 0;
		}
	}
	// End of initialization

	for (i = design->mode1; i < N; i++) // Ignore the pseudo node (pseudo node used to simulate the self-activation)
	{
		if (design->mode == 2 && i == design->mode1) // Ignore the trigger node
		{
			design->NS2[i] = 1;
			design->solution[i].curr = 1;
			if (design->solution[i].max < 1)
			{
				design->solution[i].max = 1;
				design->solution[i].edge = realloc (design->solution[i].edge, sizeof (int) * N);
			}

			for (j = 0; j < N; j++)
			{
				if (i == j)
					design->solution[i].edge[j] = -1;
				else
					design->solution[i].edge[j] = 0;
				/*
				if (design->pathway[j] == design->pathway[N + j])
					design->solution[i].edge[j] = 0;
				else if (design->pathway[j] == 1)
					design->solution[i].edge[j] = -1;
				else
					design->solution[i].edge[j] = 1;
				*/
			}
			continue;
		}

		design->min_link[i] = N;

		if (check_satisfy (design, i) == 0)
			return 0;

		design->NS *= design->NS2[i];
		design->min_links += design->min_link[i];

		if ((i == design->mode1 + 1 && design->mode == 2) || (i == design->mode1 && design->mode != 2) || design->NS2[i] < design->NS2_min)
			design->NS2_min = design->NS2[i];

		// calculate the distribution of number of links
		if (design->bdflag)
		{
			if (design->mode == 2)
				combine_dist (N, N, design->links_dist + i * (N + 1), design->Links_dist, i - design->mode1 - 1);
			else
				combine_dist (N, N, design->links_dist + i * (N + 1), design->Links_dist, i - design->mode1);
		}
        
		// calculate the degree distribution
		if (design->dgflag)
		{
			if (design->mode == 2)
				combine_dist (N, 2 * N, design->degree_dist + i * (2 * N + 1), design->Degree_dist, i - design->mode1 - 1);
			else
				combine_dist (N, 2 * N, design->degree_dist + i * (2 * N + 1), design->Degree_dist, i - design->mode1);
		}
	}

	return 1;
}
